ptm = proc.time()

+ + =

Pacotes e Leitura

Alguns pacotes R utilizados nessa rotina

library(tidyverse)          # Manipulação de banco de dados e análise exploratória
library(highcharter)        # visualização Gráfica
library(DT)                 # Construção de Tabelas 
library(knitr)              # Opções Rmarkdown (inclusão de tabelas, imagens, etc.)
library(kableExtra)         # Construção de Tabelas 
library(rpart.plot)         # Recursive Partitioning and Regression Trees
library(randomForest)       # Modelo Random Forest
library (e1071)             # Modelo Support Vector Machine
library(corrplot)           # Matriz de correlação - Visualização gráfica
library(wordcloud2)         # Numvem de palavras (Word Cloud)
library(caret)              # Matriz de confusão (Confusion Matrix) 
library(rpart)
#devtools::install_github("ropenscilabs/icon")
library(icon)

1 Introdução

Frequentemente utilizamos serviços de classificação, seja com o intuito de recomendar/sugerir algo a um usuário final, como Netflix, Mercado Livre, OLX dentre outros, em testes de diagnóstico de doenças, por exmeplo, se um paciente pode-se utilizar ummodelo de classificação como Random Forest (RF) ou Support Vector Machine (SVM) a fim de classificar pacientes como doente, não doente ou suspeito. Além de outros vários exemplos. O conjunto de dados explanados a seguir foram extraídos da plataforma Kaggle e são referentes ao catálogo de filmes da plataforma Internet Movie DataBases - IMDB, o intuito aqui é criar um sistema de classificação de estrelas (rating/score) do catálogo de avaliação de filmes do IMDB, segmentando os filmes em duas categorias de filmes: alto e baixo.

Esse projeto irá explorar o conjunto de dados a fim de gerar insights, e aplicar dois algoritmos de classificação e, posteriormente, compará-los a fim de identificar o melhor modelo que se ajusta aos dados. São eles os modelos de Random Forest - RF (baseado na clusterização de Breiman) e Support Vector Machine - SVM para classificação de filmes nas seguintes faixas/categorias de estrelas (Rating):

Categoria B: ‘Baixo índice de aprovação (low score)’: [0, 6.5) ★
Categoria A: ‘Alto índice de aprovação (high score)’: [6.5, 10) ★

Para este, se faz necessário um largo conjunto de dados que, por sua vez, deve ser bem estruturado e consistente. Portanto, é fundamental manipular, limpar e remover informações faltantes dos dados para a implementação do modelo de Machine Learning do Random Forest que será treinado com orientação a esses dados.

Além disso, esse trabalho irá levantar quais os fatores mais importantes para que um filme tenha uma alta classificação (categoria de faixas de estrelas), via modelagem de Random Forest, de tal forma que essas mesmas variáveis/fatores destacados como importantes no algoritmo do RF possam ser aproveitados posteriormente em análises futuras de classificação como no caso do classificador via SVM que será, também, aplicado aqui. Resultados e script foram gerados utilizando recursos e linguagens em \(R\), \(\LaTeX\) e \(Markdown\).

2 Banco de Dados

O conjunto de dados utilizado foi o de título IMDB 5000 extraído da plataforma kaggle. As informações contidas no banco foram catalogadas de filmes publicados ao longo de 100 anos em 66 países (entre 1916 e 2016) da plataforma IMDB - Internet Movie DataBases, o Arquivo original contém 5044 filmes (observações) e 28 variáveis descritas a seguir.

RECOMEND.METADATA = readxl::read_xlsx(path = ".../DB/IMDB.xlsx", sheet = 1)
db = RECOMEND.METADATA

Os dados em questão são públicos e disponíveis para download clicando AQUI.

3 Análise Exploratória dos dados

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

3.1 Estrutura dos dados

## O banco de dados possui 5043 observações e 28 variáveis

Excluindo algumas variáveis que não serão utilizadas.

1. As variáveis color, director_name, language, country, content_rating, em especial, foram removidas uma vez que estas não se mostraram importantes no algoritmo de importância do Random Forest de Breiman.
1. Já as variáveis director_facebook_likes, actor_3_facebook_likes, actor_2_name, actor_1_facebook_likes, actor_1_name, actor_3_name,-facenumber_in_poster, actor_2_facebook_likes, movie_imdb_link foram extraídas de forma determinística por escolha dos autores deste.

O motivo de remover as variáveis em (1) é devido à tentativa de reduzir o número de missing na base de dados. no passo 3.1 de Tratamento de NA’s.

#db=RECOMMEND.METADATA
db = db %>% select(- director_facebook_likes, -actor_3_facebook_likes, -actor_2_name, -actor_1_facebook_likes, -actor_1_name, -actor_3_name,-facenumber_in_poster, -actor_2_facebook_likes, -color, -director_name, -language, -country, -content_rating, -movie_imdb_link)

## O banco de dados agora possui 14 variáveis. Ou seja, foram removidas  14 variáveis no passo anterior.

Estrutura dos dados

glimpse(db) #str(db)

## Observations: 5,043
## Variables: 14
## $ num_critic_for_reviews    <int> 723, 302, 602, 813, NA, 462, 392, 32...
## $ duration                  <int> 178, 169, 148, 164, NA, 132, 156, 10...
## $ gross                     <int> 760505847, 309404152, 200074175, 448...
## $ genres                    <fct> Action|Adventure|Fantasy|Sci-Fi, Act...
## $ movie_title               <fct> Avatar , Pirates of the Caribbean: A...
## $ num_voted_users           <int> 886204, 471220, 275868, 1144337, 8, ...
## $ cast_total_facebook_likes <int> 4834, 48350, 11700, 106759, 143, 187...
## $ plot_keywords             <fct> avatar|future|marine|native|parapleg...
## $ num_user_for_reviews      <int> 3054, 1238, 994, 2701, NA, 738, 1902...
## $ budget                    <dbl> 237000000, 300000000, 245000000, 250...
## $ title_year                <int> 2009, 2007, 2015, 2012, NA, 2012, 20...
## $ imdb_score                <fct> 7.9, 7.1, 6.8, 8.5, 7.1, 6.6, 6.2, 7...
## $ aspect_ratio              <fct> 1.78, 2.35, 2.35, 2.35, , 2.35, 2.35...
## $ movie_facebook_likes      <int> 33000, 0, 85000, 164000, 0, 24000, 0...

3.2 Transformação de variáveis

Transformando as variáveis imdb_score e aspect_ratio em numéricas

db$imdb_score = as.numeric(as.character(db$imdb_score))
db$aspect_ratio = as.numeric(as.character(db$aspect_ratio))

hc1 = hchart(db$imdb_score, color = "#e8bb0b", name = "imdb_score") %>% 
        hc_title(text = "Histograma dos votos na plataforma IMDB") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc1

hc2 = hchart(db$num_voted_users, color = "#786eea", name = "imdb_num_votos") %>% 
        hc_title(text = "Histograma dos número de votos na plataforma IMDB") %>%
        hc_exporting(enabled = TRUE, filename = "Fig2-Pimenta"); hc2

Criando uma função para cálculo da moda

Moda <- function(x) {
     ux <- unique(x)
     ux[which.max(tabulate(match(x, ux)))]
}

## Média, moda e mediana da variável imdb_score são, respectivamente, 6.44 6.7 6.6

Extraindo valores da variável gênero e transformando em dummies

gg <- as.character(db$genres); gg <- gsub("-", "_", as.character(gg))

#t <- unlist(strsplit(gg[1],split = "\\|"))
tem1 <- data.frame() 
for(i in 1:length(gg)){
        tem <- tem1
        t <- unlist(strsplit(gg[i],split = "\\|"))
        temp <- data.frame(t)
        tem1 <- rbind(tem,temp)
}

Gen <- unique(tem1); Gen <- gsub("-", "_", as.character(Gen$t))
cat("Existem", length(Gen), "valores de gêneros únicos de filme no banco de dados.")

## Existem 26 valores de gêneros únicos de filme no banco de dados.

Genname <- Gen

fe <- matrix(data = 0, nrow = length(gg), ncol = length(Genname))
fe <- data.frame(fe); colnames(fe) <- Genname

i=j=0
for(i in 1:length(gg)){
        for(j in 1:length(Genname)){
                g <- grepl(Genname[j], gg[i])
                if(g == TRUE){
                        fe[i, j] <- 1        
                }
        }
}

NumGen = as_tibble(rbind(apply(fe,2,sum)))
NumGen = gather(NumGen, key = "variables", value = "num_gender")
NumGen = NumGen[order(NumGen$num_gender, decreasing = TRUE), ]


hc3 <- highchart() %>%
  hc_add_series(data = NumGen$num_gender, 
                type = "bar",
                name = "# de filmes",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 0, valuePrefix = "", valueSuffix = ""), color="blue") %>%
  hc_yAxis(title = list(text = "Quantitativo de filmes"), 
           allowDecimals = TRUE, max = (max(NumGen$num_gender)+103),
           labels = list(format = "{value}")) %>%
  hc_xAxis(title = list(text = "Gênero de filme"),
           categories = NumGen$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Quantitativo de filmes por gênero",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "{point.y} filmes")%>%
                 #pointFormat = "Variável: {point.x} <br> Missing: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "F3-filmes-genero-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc3; remove(gg,t,tem, temp, tem1,Gen, g, i, j)

Algumas variáveis de gêneros possuem pouquíssimos filmes serão removidas do banco. São elas Film_Noir, Short, News, Reality_TV e Game_Show.

fe$Film_Noir <- NULL
fe$Short <- NULL
fe$News <- NULL
fe$Reality_TV <- NULL
fe$Game_Show <- NULL

Unificando a base de dados com as novas variáveis de gêneros únicos descobertos.

db1 = as_tibble(cbind(as.data.frame(db),fe))
cat("Foram inseridas ", dim(fe)[2],"novas variáveis provenientes dos gêneros únicos descobertos no passo anterior.")

## Foram inseridas  21 novas variáveis provenientes dos gêneros únicos descobertos no passo anterior.

Semelhante à variável de gênero, foi feito o split por palavra chave da variável keyword, esse, por sua vez, continha 6780 palavras únicas. E então, geramos uma nuvem de palavras chave, com o auxílio do pacote wordcloud2 com o seguinte comando, após obter o data frame do nome das palavras chaves e frequência.

library(wordcloud2)
wordcloud2(NumKeyWord)

3.3 Filtrando filmes a partir de 1980

Excluindo filmes lançados antes de 1980

hc00 = hchart(db$title_year, color = "#53a074", name = "imdb1_score") %>% 
        hc_title(text = "Histograma do ano de publicação dos filmes") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc00

Percebemos pelo gráfico anterior que existem poucos filmes publicados antes de 1980 e estes podem não ser representativos. Decidimos, então, por trabalhar apenas com os filmes que foram publicados a aprtir de 1980

db <- db[db$title_year >= 1980,]
hc01 = hchart(db$title_year, color = "#79d8a2", name = "imdb1_score") %>% 
        hc_title(text = "Histograma do ano de publicação dos filmes") %>%
        hc_exporting(enabled = TRUE, filename = "Fig1-Pimenta"); hc01

3.4 Tratamento de NA, zeros e duplicatas

Porcentagem de NA por variável

db_miss <- db %>% summarise_all(funs(sum(is.na(.))/n()))
db_miss <- gather(db_miss, key = "variables", value = "percent_missing")
db_miss$percent_missing = 100*db_miss$percent_missing
db_miss = db_miss[order(db_miss$percent_missing, decreasing = TRUE), ]
#db_miss

hc4 <- highchart() %>%
  hc_add_series(data = db_miss$percent_missing, 
                type = "bar",
                name = "Porcentagem de missing",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = " %")) %>%
  hc_yAxis(title = list(text = "Porcentagem de missing"), 
           allowDecimals = TRUE, max = 20,
           labels = list(format = "{value}%")) %>%
  hc_xAxis(categories = db_miss$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Porcentagem de missinig por variável",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Missing: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Missing: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "Fig0-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc4

Removendo todas as observações que contêm NA.

db1 <- db1 %>% drop_na(gross, budget, aspect_ratio 
,title_year 
,num_critic_for_reviews, num_user_for_reviews)
#glimpse(db1)

db1 <- db1 %>% na.omit()
db1 <- db1 %>% drop_na()

Removendo dados duplicados

cat("Existem", sum(duplicated(db1)), "filmes duplicados na base de dados.")

## Existem 33 filmes duplicados na base de dados.

db1 <- db1[!duplicated(db1), ]

## O novo banco de dados, sem observações faltantes, possui 3782 observações. Ou seja, no processo de remoção de valores faltantes foram perdidas 1000 observações. Finalmente, removemos todas as observações duplicadas e faltantes do banco.

Porcentagem de ZEROS por variável

zeros <- (colSums(db1==0)/nrow(db1)*100); var <- names(zeros)
db_zero <- data.frame(var,zeros); rownames(db_zero) <- NULL
db_zero <- db_zero[order(db_zero$zeros, decreasing = TRUE), ]

hc4_1 <- highchart() %>%
  hc_add_series(data = db_zero$zeros, 
                type = "bar",
                name = "Porcentagem de zeros",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = " %"), color="pink") %>%
  hc_yAxis(title = list(text = "Porcentagem de zero"), 
           allowDecimals = TRUE, max = 100,
           labels = list(format = "{value}%")) %>%
  hc_xAxis(categories = db_zero$var,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Porcentagem de zeros por variável",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Zeros: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Missing: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "Fig00-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc4_1

Aqui podemos perceber que os zeros estão concentrados nas variáveis de gênero que são binárias e não faz sentido analizar os zeros nessas variáveis.

paste0("Entretanto, pode-se perceber uma proporção de zeros na variável movie_facebook_likes igual a ",round(db_zero[22,2],2), "% e, também na variável cast_total_facebook_likes de ",round(db_zero[23,2],2),"%. Decidimos por remover a variável movie_facebook_likes.")

## [1] "Entretanto, pode-se perceber uma proporção de zeros na variável movie_facebook_likes igual a 46.17% e, também na variável cast_total_facebook_likes de 0.13%. Decidimos por remover a variável movie_facebook_likes."

db1$movie_facebook_likes <- NULL

3.5 Weighted Rating (WR) - IMDB_score

Um passo importante é penalizar a variável de escore imdb_score pelo número de votos recebidos. Ver mais no estimador de Shrinkage

\(WR = \frac{v}{v+m} \times R + \frac{m}{v+m} \times C\)

Onde,

R = as.numeric(db1$imdb_score)
v = as.numeric(db1$num_voted_users)
m = summary(db1$num_voted_users)[2] # 1st Qu.
C = mean(db1$imdb_score)
db1$WR = (v/(v+m))*R + (m/(v+m))*C

R = Escore médio dos votos para o título do filme dado pelos usuários do IMDB = (imdb_score)
v = Número de usuários que votaram = (num_voted_users)
m = Mínimo de votos requerido (atualmente 7.000)
C = O escore médio de todos os 3766 filmes (atualmente 6,5)

Medidas pontuais e de dispersão de imdb_score e WR

summary(db1$imdb_score)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.600   5.900   6.600   6.468   7.200   9.300

summary(db1$WR)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.578   6.166   6.498   6.573   6.983   9.269

cat("Além disso, conseguimos reduzir o desvio padrão que orignalmente era de ",round(sd(db1$imdb_score),2),",para a variável imdb_score, e agora passou a ser ",round(sd(db1$WR),2), "para a variável WR.")

## Além disso, conseguimos reduzir o desvio padrão que orignalmente era de  1.05 ,para a variável imdb_score, e agora passou a ser  0.72 para a variável WR.

Selecionando uma amostra aleatória (a.a.) de tamanho \(n=800\) e representando em um gráfico de dispersão com valores reais vs ajustados.

set.seed(123654)
amostra0 = sample(x = 1:dim(db1)[1], size = 800, replace = FALSE)
dbX = db1[amostra0,] %>% 
  select(movie_title,num_voted_users,imdb_score, WR)
  
dss <- map(c("cross"), function(s){
  
  x <- as.numeric(dbX$imdb_score)
  y <- as.numeric(dbX$WR)
  
  list(name = s,
       data = list_parse(data_frame(x, y)),
       marker = list(symbol = s, enabled = TRUE), lineColor = "#56667a")
  
})
#dss[[1]]$data[amostra1]

hc5 = highchart() %>% 
  hc_chart(type = "scatter", color = "#56667a") %>% 
  hc_title(text = "Score IMDB vs WR (calculado pelo estimador de Shrinkage)") %>%
  hc_subtitle(text = "800 filmes selecionados via amostra aleatória simples") %>%
  hc_xAxis(title = list(text = "x: imdb_score"), 
           allowDecimals = TRUE, labels = list(format = "{value}★")) %>%
  hc_yAxis(title = list(text = "y: WR (calibrado)"),
           allowDecimals = TRUE, labels = list(format = "{value}★")) %>%
  hc_exporting(enabled = TRUE, filename = "F3-Pimenta") %>%
  hc_add_series_list(dss); hc5; remove(dbX, dss)

Análise do cálculo de IMDB

Para filmes com # de votos recebidos MENOR que 18mil (m: votos mínimos requeridos):
- imdb_score > 6,5: DECRESCIMENTO
- imdb_score < 6,5: CRESCIMENTO
Para filmes com # de votos recebidos MAIOR que 18mil (m: votos mínimos requeridos):
- imdb_score > 6,5: DECRESCIMENTO
- imdb_score < 6,5: CRESCIMENTO

Ou seja, para os filmes catalogados com m muito inferior a 18mil o novo escore calibrado teve maior diferentça que aqueles superior a 18 mil. Além disso, quando scores são maiores que 6,5 o WR tende a cair, caso contrário o valor pode aumentar. Quanto maior o número de votos recebidos, menor a diferença do valor de IMDB_score e WR.

set.seed(123654)
amostra2 = sample(x = 1:dim(db1)[1], size = 35, replace = FALSE)
db1[amostra2,] %>% 
  select(movie_title,num_voted_users,imdb_score, WR, budget) %>% 
  datatable(rownames = amostra2,options = list(searchin = TRUE, scrollX = TRUE, pageLength = 5))

3.6 Visualizações Finais

Finalmente, segue uma a.a.da base de dados após limpeza e análise exploratória.

set.seed(123851)
amostra3 = sample(x = 1:dim(db1)[1], size = 35, replace = FALSE)
db1[amostra3,] %>%
datatable(rownames = amostra3,options = list(searchin = TRUE, scrollX = TRUE, pageLength = 5))

Para os dados quantitativos a seguinte Matriz de correlação. Foi gerada utilizando o pacote corrplot

#install.packages("corrplot")
library(corrplot)  
res1 <- cor.mtest(db_cor, conf.level = .95)
res2 <- cor.mtest(db_cor, conf.level = .99)

corrplot::corrplot(cor(db_cor), method = "color", col = col(200),
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90, # Text label color and rotation
         # Combine with significance
         p.mat = res1$p, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag = FALSE)

4 Preparação e Partição dos Dados

4.1 Preparação de dados

dist1 = 2.5
db2 = db1 %>% filter(WR >= (mean(db1$WR)- dist1) & WR <= (mean(db1$WR)+ dist1))
paste('Excluímos ', dim(db1)[1] - dim(db2)[1], ' filmes da base. O critério aqui utilizado foi o de permanecer, apenas, filmes distantes em até ', dist1, ' estrelas da média de estrelas ',round(mean(db1$WR),2),'.')

## [1] "Excluímos  16  filmes da base. O critério aqui utilizado foi o de permanecer, apenas, filmes distantes em até  2.5  estrelas da média de estrelas  6.57 ."

par(mfrow=c(2,2))
plot(db2$imdb_score, ylim = c(0,10), main = "Plot da variável imdb_score")
abline(h=mean(db2$imdb_score), col = "purple")
plot(db2$WR, ylim = c(0,10), main = "Plot da variável WR (imdb_score penalizado)")
abline(h=mean(db2$WR), col = "purple")
hist(db2$imdb_score, main = "Histograma da variável imdb_score", probability = TRUE, xlab = "Estrelas", ylab = "Porcentagem de Estrelas")
#curve(dnorm(x, mean=mean(db2$imdb_score), sd=sd(db2$imdb_score)), add=TRUE, col = "red")
points(seq(min(db2$imdb_score), max(db2$imdb_score), length.out = length(db2$imdb_score)), 
           dnorm(seq(min(db2$imdb_score), max(db2$imdb_score), length.out = length(db2$imdb_score)), 
                     mean(db2$imdb_score), sd(db2$imdb_score)), type = "l", col = "red")
lines(density(db2$imdb_score), col="blue", lwd=1) # add a density estimate with defaults
lines(density(db2$imdb_score, adjust=2), lty="dotted", col="darkgreen", lwd=1) 
hist(db2$WR, main = "Histograma da variável WR (imdb_score penalizado)", probability = TRUE, xlab = "Estrelas", ylab = "Porcentagem de Estrelas", ylim = c(0,.8))
points(seq(min(db2$WR), max(db2$WR), length.out = length(db2$WR)), 
           dnorm(seq(min(db2$WR), max(db2$WR), length.out = length(db2$WR)), 
                     mean(db2$WR), sd(db2$WR)), type = "l", col = "red")
lines(density(db2$WR, ), col="blue", lwd=1) # add a density estimate with defaults
lines(density(db2$WR, adjust=2), lty="dotted", col="darkgreen", lwd=1)

Transformação da variável contínua WRem categórica WR_Grp

Transformando a variável WR em fator, com as seguintes categorias:

Categoria ‘Baixo índice de aprovação (low score)’: [0, 6.5) ★
Categoria ‘Alto índice de aprovação (high score)’: [6.5, 10) ★

movie = db2
Grp <- function(tn){
  tn = abs(tn)
    if (tn >= 0 & tn < median(db2$WR)){
        return(paste0('[0, ',round(median(db2$WR),2),')'))
    }else if(tn >= median(db2$WR)){
        return(paste0('[',round(median(db2$WR),2),', 10]'))
    }
}
# apply the Group function to the WR column
movie$WR_Grp <- sapply(movie$WR,Grp)
# set as factor the new column
movie$WR_Grp <- as.factor(movie$WR_Grp)
#View(head(movie))
table(movie$WR_Grp)

## 
##  [0, 6.5) [6.5, 10] 
##      1883      1883

# apply the Group function to the WR column
imdb_score_Grp <- sapply(movie$imdb_score,Grp)
# set as factor the new column
imdb_score_Grp <- as.factor(imdb_score_Grp)
#View(head(movie))
table(imdb_score_Grp)

## imdb_score_Grp
##  [0, 6.5) [6.5, 10] 
##      1696      2070

Como podemos ver, na base sem a transformação e sem a exclusão de filmes com mais de 2,5 estrelas distantes da média tinha quase 1700 filmes com baixa classificação contra 2070 com alta classificação.

O critério de selecionar filmes distantes em até 2,5 estrelas da média de 6,57 juntamente com a transoformação da variável pelo cálculo de WR permitiu construir uma nova base de dados com categorias (Alta e Baixa) muito mais balanceadas, com um número de filmes em ambas as categorias igual a 1883.

Remoção de variáveis Remove as variáveis imdb_score, genres, plot_keywords, movie_imdb_link e WR.

movie$imdb_score <- NULL
movie$genres <- NULL
movie$plot_keywords <- NULL
movie$movie_imdb_link <- NULL
movie$WR <- NULL

#str(as_tibble(movie))
glimpse(as_tibble(movie))

## Observations: 3,766
## Variables: 32
## $ num_critic_for_reviews    <int> 723, 302, 602, 813, 462, 392, 324, 6...
## $ duration                  <int> 178, 169, 148, 164, 132, 156, 100, 1...
## $ gross                     <int> 760505847, 309404152, 200074175, 448...
## $ movie_title               <fct> Avatar , Pirates of the Caribbean: A...
## $ num_voted_users           <int> 886204, 471220, 275868, 1144337, 212...
## $ cast_total_facebook_likes <int> 4834, 48350, 11700, 106759, 1873, 46...
## $ num_user_for_reviews      <int> 3054, 1238, 994, 2701, 738, 1902, 38...
## $ budget                    <dbl> 237000000, 300000000, 245000000, 250...
## $ title_year                <int> 2009, 2007, 2015, 2012, 2012, 2007, ...
## $ aspect_ratio              <dbl> 1.78, 2.35, 2.35, 2.35, 2.35, 2.35, ...
## $ Action                    <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, ...
## $ Adventure                 <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ Fantasy                   <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, ...
## $ Sci_Fi                    <dbl> 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, ...
## $ Thriller                  <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Documentary               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Romance                   <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, ...
## $ Animation                 <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Comedy                    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Family                    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, ...
## $ Musical                   <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Mystery                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
## $ Western                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Drama                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ History                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Sport                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Crime                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Horror                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ War                       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Biography                 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Music                     <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ WR_Grp                    <fct> [6.5, 10], [6.5, 10], [6.5, 10], [6....

cat("A base de dados que iremos trabalhar no modelo tem", dim(movie)[1] , "observações e ", dim(movie)[2]-1, "variáveis, sem contar a variável que contém os títulos dos filmes.")

## A base de dados que iremos trabalhar no modelo tem 3766 observações e  31 variáveis, sem contar a variável que contém os títulos dos filmes.

4.2 Partição dos dados

Particionando a base de dados em Treino e Teste, esses dois (Treino e Teste) também terão armazenos os nomes dos filmes selecionados via amostragem probabilística dos dados originais separadamente das bases de Treino e Teste.

Posteriormente, removemos os rótulos dos filmes nas bases Treino e Teste

set.seed(9182345)
ind <- sample(2, nrow(movie), replace = T, prob = c(0.7, 0.3))
train <- movie[ind==1, -4]
test <- movie[ind==2, -4]
trainMovie <- movie[ind==1, 4]
testMovie <- movie[ind==2, 4]

Porcentagem da Distribuição dos filmes por categorias nas bases de treino e teste

round(table(train$WR_Grp)/sum(table(train$WR_Grp))*100,2)

## 
##  [0, 6.5) [6.5, 10] 
##     49.79     50.21

round(table(test$WR_Grp)/sum(table(test$WR_Grp))*100,2)

## 
##  [0, 6.5) [6.5, 10] 
##      50.5      49.5

hchart(train$WR_Grp, colorByPoint = TRUE, name = "Escore - Base Treino")

hchart(test$WR_Grp, colorByPoint = TRUE, name = "Escore - Base Teste")

MODELAGEM

5 Modelagem 1 - Random Forest (RF)

5.1 Random Forest (RF) - Metodologia

Descrição 1. Random Forest foi desenvolvido para agregar árvores de decisão (modelo de classificação);
2. Pode ser usado para modelo de classificação (p/ var. resposta categórica) ou regressão (no caso de haver variável resposta contínua);
3. Evita overfitting;
4. Permite trabalhar com um largo número de características de um conjunto de dados;
5. Auxilia na seleção de variáveis baseada em um algoritmo que calcula a importância por variável (assim, tendo conhecimento de quais variáveis são mais importantes, podemos usar essa informação para outros modelos de classificação);
6. User-friendly: apenas 2 parâmetros livres:

Trees - ntrees, default 500 (Nº de árvores);
Variáveis selecionadas via amostragem aleatória candidatas à cada “split” (quebra da árvore) - mtry, default \(\sqrt{p}\) p/ classificação e \(\frac{p}{3}\) p/ regressão (p: nº de features/variáveis);

Passo-a-Passo

É realizado em 3 passos:

Desenha as amostras via bootstrap do número de árvores ntrees;
Para cada amostra via bootstrap, cresce o número de árvores “un-puned” para a escolha da melhor quebra da árvore baseado na amostra aleatória do valor predito de mtry a cada nó da árvore;

1. Faz classificação de novos valores usando a maioria de votos p/ classificação e usa a média p/ regressão baseada nas amostras de ntrees.

Exemplo

5.2 Random Forest - Aplicação e Resultados

Inicialmente utilizaremos o pacote randomForest que implmenta o algoritmo de Random Forest de Breiman (baseado na clusterização de Breiman, originalmente codificada em Fortran) que tem por finalidade classificar e/ou criar regressão. Além disso, pode ser usado em um modelo não supervisionado para avaliar proximidades entre pontos.

Estamos usando, a partir daqui, a base de treino.

#library(randomForest)
#library(rpart)
#library(rpart.plot)
#rf <- randomForest(proximity = T,ntree = 38,do.trace = T,WR~.,data=training)
set.seed(9984512)
rf <- randomForest(WR_Grp~num_voted_users,data=train)
rf

## 
## Call:
##  randomForest(formula = WR_Grp ~ num_voted_users, data = train) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 40.42%
## Confusion matrix:
##           [0, 6.5) [6.5, 10] class.error
## [0, 6.5)       788       540   0.4066265
## [6.5, 10]      538       801   0.4017924

attributes(rf)

## $names
##  [1] "call"            "type"            "predicted"      
##  [4] "err.rate"        "confusion"       "votes"          
##  [7] "oob.times"       "classes"         "importance"     
## [10] "importanceSD"    "localImportance" "proximity"      
## [13] "ntree"           "mtry"            "forest"         
## [16] "y"               "test"            "inbag"          
## [19] "terms"          
## 
## $class
## [1] "randomForest.formula" "randomForest"

Olhando as 6 primeiras observações real X predito

p1 <- predict(rf,train)
head(p1)

##         1         2         3         4         5         6 
## [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10]  [0, 6.5) 
## Levels: [0, 6.5) [6.5, 10]

head(train$WR_Grp)

## [1] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [0, 6.5) 
## Levels: [0, 6.5) [6.5, 10]

Matriz de confusão

library(caret)
library(e1071)
confusionMatrix(p1, train$WR_Grp)

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  [0, 6.5) [6.5, 10]
##   [0, 6.5)      1325        13
##   [6.5, 10]        3      1326
##                                           
##                Accuracy : 0.994           
##                  95% CI : (0.9903, 0.9966)
##     No Information Rate : 0.5021          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.988           
##  Mcnemar's Test P-Value : 0.02445         
##                                           
##             Sensitivity : 0.9977          
##             Specificity : 0.9903          
##          Pos Pred Value : 0.9903          
##          Neg Pred Value : 0.9977          
##              Prevalence : 0.4979          
##          Detection Rate : 0.4968          
##    Detection Prevalence : 0.5017          
##       Balanced Accuracy : 0.9940          
##                                           
##        'Positive' Class : [0, 6.5)        
##

RF_importance = randomForest::importance(rf)[order(randomForest::importance(rf)[,1], decreasing = TRUE), ]
knitr::kable(RF_importance)

x
1323.525

Verificamos que as variávei Game_Show, Sci_Fi, Reality_TV, News e Film_Noir não foram relevantes para o algoritmo do random forest.

randomForest::varImpPlot(rf)

Taxa de Erro - Random Forest

plot(rf)
legend('topright', colnames(rf$err.rate), col=1:5, fill=1:5)

Observamos que, a partir do número de árvores geradas \(ntrees > 300\) o erro OOB (Out of Bag) não pode ser melhorado.

Ajustando e melhorando estimativas

Além disso, iremos alterar alguns parâmetros da função randomForest como o número de ntrees e mtry. Assim, repetimos o algoritmo do Random Forest, ainda usando a base treino.

# Tune mtry
x = as.data.frame(train1[,-31])
y = (as.factor(train1$WR_Grp))
t <- tuneRF(x = x, y = y,
       stepFactor = 0.3,
       plot = TRUE,
       ntreeTry = 300,
       trace = TRUE,
       improve = 0.05)

## mtry = 5  OOB error = 18.64% 
## Searching left ...
## mtry = 17    OOB error = 18.52% 
## 0.006036217 0.05 
## Searching right ...
## mtry = 1     OOB error = 23.81% 
## -0.277666 0.05

Aparentemente, 5 é um bom candidato ao valor de \(m_{try}\)

#set.seed(093180)
set.seed(998451)
rf1 <- randomForest(WR_Grp~.,data=train1, 
                    ntree = 300, 
                    mtry = 5, 
                    importance = TRUE,
                    proximity = TRUE)
rf1

## 
## Call:
##  randomForest(formula = WR_Grp ~ ., data = train1, ntree = 300,      mtry = 5, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 300
## No. of variables tried at each split: 5
## 
##         OOB estimate of  error rate: 18.45%
## Confusion matrix:
##           [0, 6.5) [6.5, 10] class.error
## [0, 6.5)      1105       223   0.1679217
## [6.5, 10]      269      1070   0.2008962

RF_importance1 = randomForest::importance(rf1)[order(randomForest::importance(rf1)[,1], decreasing = TRUE), ]
plot(rf1)
legend('topright', colnames(rf$err.rate), col=1:5, fill=1:5)

Removendo variáveis com pouca importância

Decidimos por retirar todos os fatores que retornaram importância abaixo de 10 em MeanDecreaseGini. Portanto, ficaremos apenas com as seguintes variáveis:

RF1 = data.frame(variables = rownames(RF_importance1), importance = RF_importance1[,4])
RF1 = RF1[order(RF1$importance, decreasing = TRUE),]
rownames(RF1) <- NULL
RF1 = RF1[1:10,]
hc6_1 <- highchart() %>%
  hc_add_series(data = RF1$importance, 
                type = "bar",
                name = "Importância",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = "")) %>%
  hc_yAxis(title = list(text = "Importância"), 
           allowDecimals = TRUE, max = 200,
           labels = list(format = "{value}")) %>%
  hc_xAxis(title = list(text = "Fatores"),
           categories = RF1$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Importância por fator - Random Forest",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Importância: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Importância: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "F6_1-importance-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc6_1

Selecionando variáveis importantes

train2 = train %>% select(num_voted_users, gross, duration, num_user_for_reviews, budget, cast_total_facebook_likes, title_year, num_critic_for_reviews, Drama, Horror, aspect_ratio, Action, Comedy, WR_Grp) %>% droplevels()
test2 = test %>% select(num_voted_users, gross, duration, num_user_for_reviews, budget, cast_total_facebook_likes, title_year, num_critic_for_reviews, Drama, Horror, aspect_ratio, Action, Comedy, WR_Grp) %>% droplevels()
movie = rbind(train2, test2)

# Tune mtry
x = as.data.frame(train2[,-14])
y = (as.factor(train2$WR_Grp))
t <- tuneRF(x = x, y = y,
       stepFactor = .7,
       plot = TRUE,
       ntreeTry = 150,
       trace = TRUE,
       improve = 0.05)

## mtry = 3  OOB error = 20.25% 
## Searching left ...
## mtry = 5     OOB error = 20.32% 
## -0.003703704 0.05 
## Searching right ...
## mtry = 2     OOB error = 20.85% 
## -0.02962963 0.05

Modelo Final Os parâmetros utilizados, finalmente serão \(m_{try} = 3\) e \(n_{tree} = 150\)

set.seed(093180)
rf_final <- randomForest(WR_Grp~.,data = train2,
                    ntree = 150, 
                    mtry = 3, 
                    importance = TRUE,
                    proximity = TRUE)
rf_final; #attributes(rf_final)

## 
## Call:
##  randomForest(formula = WR_Grp ~ ., data = train2, ntree = 150,      mtry = 3, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 150
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 20.02%
## Confusion matrix:
##           [0, 6.5) [6.5, 10] class.error
## [0, 6.5)      1089       239   0.1799699
## [6.5, 10]      295      1044   0.2203137

RF_importance1 = randomForest::importance(rf_final)[order(randomForest::importance(rf_final)[,1], decreasing = TRUE), ]
plot(rf_final)
legend('topright', colnames(rf$err.rate), col=1:5, fill=1:5)

Predição e matriz de confusão - train data

library(caret)
p1 <- predict(rf_final,train2)
head(p1)

##         1         2         3         4         5         6 
## [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10]  [0, 6.5) 
## Levels: [0, 6.5) [6.5, 10]

head(train2$WR_Grp)

## [1] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [0, 6.5) 
## Levels: [0, 6.5) [6.5, 10]

confusionMatrix(p1, train2$WR_Grp)

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  [0, 6.5) [6.5, 10]
##   [0, 6.5)      1328         0
##   [6.5, 10]        0      1339
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9986, 1)
##     No Information Rate : 0.5021     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##  Mcnemar's Test P-Value : NA         
##                                      
##             Sensitivity : 1.0000     
##             Specificity : 1.0000     
##          Pos Pred Value : 1.0000     
##          Neg Pred Value : 1.0000     
##              Prevalence : 0.4979     
##          Detection Rate : 0.4979     
##    Detection Prevalence : 0.4979     
##       Balanced Accuracy : 1.0000     
##                                      
##        'Positive' Class : [0, 6.5)   
##

Predição e matriz de confusão - test data

p2 <- predict(rf_final,test2)
head(p2)

##         1         2         3         4         5         6 
## [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10] 
## Levels: [0, 6.5) [6.5, 10]

head(test2$WR_Grp)

## [1] [6.5, 10] [0, 6.5)  [6.5, 10] [6.5, 10] [6.5, 10] [6.5, 10]
## Levels: [0, 6.5) [6.5, 10]

confusionMatrix(p2, test2$WR_Grp)

## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  [0, 6.5) [6.5, 10]
##   [0, 6.5)       469       105
##   [6.5, 10]       86       439
##                                           
##                Accuracy : 0.8262          
##                  95% CI : (0.8025, 0.8482)
##     No Information Rate : 0.505           
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6523          
##  Mcnemar's Test P-Value : 0.1928          
##                                           
##             Sensitivity : 0.8450          
##             Specificity : 0.8070          
##          Pos Pred Value : 0.8171          
##          Neg Pred Value : 0.8362          
##              Prevalence : 0.5050          
##          Detection Rate : 0.4268          
##    Detection Prevalence : 0.5223          
##       Balanced Accuracy : 0.8260          
##                                           
##        'Positive' Class : [0, 6.5)        
##

Nº de nós nas árvores

hist(treesize(rf_final),
     main = "Nº de nós por ávore",
     col = "green")

Extração de uma única árvore \(Árvore \space n_{tree}=1\)

datatable(getTree(rf_final, 1, labelVar = TRUE),  
          options = list(searchin = TRUE, pageLength = 5))

\(Árvore \space n_{tree}=149\)

datatable(getTree(rf_final, 149, labelVar = TRUE),  
          options = list(searchin = TRUE, pageLength = 5))

Gráfico de escala multidimensional da matriz de proximidade Classical multidimensional scaling (MDS) of a data matrix. Also known as principal coordinates analysis (Gower, 1966).

Note that because of numerical errors the computed eigenvalues need not all be non-negative, and even theoretically the representation could be in fewer than n - 1 dimensions.

edit(MDSplot)
fig.align="center"}
(MDIM_treino = MDSplot(rf_final, train2$WR_Grp, pch=20,  keep.forest=FALSE))v
(MDIM_teste = MDSplot(rf_final, test2$WR_Grp, pch=20,  keep.forest=FALSE))
sum(MDIM_treino$eig[1:2])
sum(MDIM_teste$eig[1:2])

Gráfico multidimensional com a variável resposta da base de TREINO.

Gráfico multidimensional com a variável resposta da base de TESTE.

Além disso, essa representação, multidimensional, explica cerca 70% da variabilidade total.

Plot do modelo rpart - Recursive Partitioning and Regression Trees - personalizando automaticamente a partir do gráfico para o tipo de resposta do modelo.

rp <- rpart::rpart(formula = WR_Grp~.,data=test2)
#rpart::plotcp(rp)
rpart.plot(rp)

rpart.plot.version1(rp)

5.3 Random Forest - Considerações

O Random Forest apresenta overfiting para os dados de treinamento, como é possível observar na matriz de confusão obtida, e a partir de 150 arvores de regressão no algoritmo, tem-se um ganho minimo ao adicionar novas árvores.

É possível comparar nossos resultados com os resultados da usuária Yueming (disponíveis no site Kaggle clicando AQUI e AQUI), sendo que a partir de seus resultados é possível obter algumas sugestões de melhora do nosso algoritmo, principalmente quanto ao tratamento dos dados (NA e zeros). Em nosso algoritmo, inicialmente, os filmes repetidos não haviam sido considerados, a autora observou 43 observações com repetição e variáveis com excesso de respostas nulas. Quanto as observações repetidas, replicamos o passo sugerido por ela e foram excluídos todos os filmes com repetição.

Em relação as respostas nulas em excesso, ela tratou essas, como NA, entretanto durante as etapas de preparação de dados, descobrimos que a única variável com tal problema era a “movie_facebook_likes” que acabou sendo excluída de nossa análise.

Ainda foram observadas poucos filmes antes de 1980. Tal situação foi enfrentada pela autora com a exclusãos dessas observações raras. Replicamos esse passo também. De uma forma geral, a autora comparou três modelos (KNN, Árvore de Decisão e Random Forest) apenas levando em consideração a acurácia obtida nos modelos.

Aqui, tentamos levar em consideração os valores, principalmente, da matriz de confusão e não somente a acurácia. E podemos concluir que:

O nosso modelo é muito bom para classificar filmes com alta e baixa aprovação. Utilizando a base de teste, o modelo retornou através da matriz de confusão uma sensibilidade de cerca de 85% e uma especificidade de 80%. Além disso, uma acurácia obtida de 83%.
É valido ressaltar, também, que o erro OOB estimado foi de 20%.
Ainda através da matriz de confusão, o nosso modelo erra em uma menor quantidade de vezes uma classificação de filmes de até 6,5★ (baixa aprovação). A estimativa de erro de classificação desses filmes é de 18% versus 22% de erro de classificação em filmes com alta aprovação.
O quadro a seguir apresenta o número de filmes em cada classificação segundo a base de teste, treino e inteira (movie). Aparentemente as categorias são bem balanceadas.

Base de dados	[0, 6.5)	[6.5, 10]
Movie (treino + teste)	1883	1883
Treino	1328	1339
Teste	555	544

Ainda conseguimos obter quais são os fatores/variáveis mais importantes para análise de classificação, dentre as variáveis existentes. Abaixo citamos as top 10 variáveis mais importantes.

RF1 = RF1[1:10,]
hc6_1 <- highchart() %>%
  hc_add_series(data = RF1$importance, 
                type = "bar",
                name = "Importância",
                showInLegend = FALSE,
                tooltip = list(valueDecimals = 2, valuePrefix = "", valueSuffix = ""),
                color="orange") %>%
  hc_yAxis(title = list(text = "Importância"), 
           allowDecimals = TRUE, max = 200,
           labels = list(format = "{value}")) %>%
  hc_xAxis(title = list(text = "Fatores"),
           categories = RF1$variables,
           tickmarkPlacement = "on",
           opposite = FALSE) %>%
  hc_title(text = "Importância por fator - Random Forest",
           style = list(fontWeight = "bold")) %>% 
  hc_subtitle(text = paste("")) %>%
      hc_tooltip(valueDecimals = 2,
                 pointFormat = "Importância: {point.y}")%>%
                 #pointFormat = "Variável: {point.x} <br> Importância: {point.y}") 
      hc_credits(enabled = TRUE, 
                 text = "Fonte: IMDB/KAGGLE. Elaboração: Ewerson Pimenta.",
                 style = list(fontSize = "10px")) %>%
  hc_exporting(enabled = TRUE, filename = "F6_1-importance-Pimenta")
#hc <- hc %>% 
#  hc_add_theme(hc_theme_darkunica())
hc6_1

6 Modelagem 2 - Support Vector Machine (SVM)

Inicialmente o modelo SVM terá como premissa as mesmas variáveis utilizadas e preparadas anteriormente no modelo Final do RF.

6.1 Ajuste do Modelo - Escolha do Kernel

Ajustando o SVM

library(e1071)
set.seed(093180)
svm_01 <- svm(WR_Grp~.,data = train2,)
#svm_01; #attributes(svm_01)
#svm_01$nSV
summary(svm_01)

## 
## Call:
## svm(formula = WR_Grp ~ ., data = train2, )
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.07692308 
## 
## Number of Support Vectors:  1528
## 
##  ( 778 750 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  [0, 6.5) [6.5, 10]

plot(svm_01, data = train2, num_voted_users~duration, 
     slice = list(title_year = 3, gross = 4))

plot(svm_01, data = train2, duration~num_voted_users, 
     slice = list(num_user_for_reviews = 3, gross = 4))

Matriz de Confusão e Erro de Classificação (Missclassification Error)

pred.svm1 = predict(svm_01, train2)
(tab.svm1 = table(Predito = pred.svm1, Real = train2$WR_Grp))

##            Real
## Predito     [0, 6.5) [6.5, 10]
##   [0, 6.5)      1119       335
##   [6.5, 10]      209      1004

(erro.kernelRAD = 1 - sum(diag(tab.svm1))/sum(tab.svm1))

## [1] 0.2039745

Ajustando SVM com kernel Linear

library(e1071)
set.seed(093180)
svm_02 <- svm(WR_Grp~.,data = train2, kernel = "linear")
summary(svm_02)

## 
## Call:
## svm(formula = WR_Grp ~ ., data = train2, kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.07692308 
## 
## Number of Support Vectors:  1481
## 
##  ( 741 740 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  [0, 6.5) [6.5, 10]

plot(svm_02, data = train2, num_voted_users~duration, 
     slice = list(title_year = 3, gross = 4))

plot(svm_02, data = train2, duration~num_voted_users, 
     slice = list(num_user_for_reviews = 3, gross = 4))

Temos que as seguintes matriz de confusão e erro de classificação para SVM com kernel Linear

pred.svm2 = predict(svm_02, train2)
(tab.svm2 = table(Predito = pred.svm2, Real = train2$WR_Grp))

##            Real
## Predito     [0, 6.5) [6.5, 10]
##   [0, 6.5)      1056       352
##   [6.5, 10]      272       987

(erro.kernelLIN = 1 - sum(diag(tab.svm2))/sum(tab.svm2))

## [1] 0.2339708

Ajustando SVM com kernel Polinomial

library(e1071)
set.seed(093180)
svm_03 <- svm(WR_Grp~.,data = train2, kernel = "polynomial")
summary(svm_03)

## 
## Call:
## svm(formula = WR_Grp ~ ., data = train2, kernel = "polynomial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  polynomial 
##        cost:  1 
##      degree:  3 
##       gamma:  0.07692308 
##      coef.0:  0 
## 
## Number of Support Vectors:  1714
## 
##  ( 856 858 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  [0, 6.5) [6.5, 10]

plot(svm_03, data = train2, num_voted_users~duration, 
     slice = list(title_year = 3, gross = 4))

plot(svm_03, data = train2, duration~num_voted_users, 
     slice = list(num_user_for_reviews = 3, gross = 4))

Temos que as seguintes matriz de confusão e erro de classificação para SVM com kernel Polynomial

pred.svm3 = predict(svm_03, train2)
(tab.svm3 = table(Predito = pred.svm3, Real = train2$WR_Grp))

##            Real
## Predito     [0, 6.5) [6.5, 10]
##   [0, 6.5)      1055       294
##   [6.5, 10]      273      1045

(erro.kernelPOLY = 1 - sum(diag(tab.svm3))/sum(tab.svm3))

## [1] 0.2125984

Ajustando SVM com kernel Sigmoidal

library(e1071)
set.seed(093180)
svm_04 <- svm(WR_Grp~.,data = train2, kernel = "sigmoid")
summary(svm_04)

## 
## Call:
## svm(formula = WR_Grp ~ ., data = train2, kernel = "sigmoid")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  sigmoid 
##        cost:  1 
##       gamma:  0.07692308 
##      coef.0:  0 
## 
## Number of Support Vectors:  1339
## 
##  ( 670 669 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  [0, 6.5) [6.5, 10]

plot(svm_04, data = train2, num_voted_users~duration, 
     slice = list(title_year = 3, gross = 4))

plot(svm_04, data = train2, duration~num_voted_users, 
     slice = list(num_user_for_reviews = 3, gross = 4))

Temos que as seguintes matriz de confusão e erro de classificação para SVM com kernel Sigmoid

pred.svm4 = predict(svm_04, train2)
(tab.svm4 = table(Predito = pred.svm4, Real = train2$WR_Grp))

##            Real
## Predito     [0, 6.5) [6.5, 10]
##   [0, 6.5)       933       426
##   [6.5, 10]      395       913

(erro.kernelSIGMOID = 1 - sum(diag(tab.svm4))/sum(tab.svm4))

## [1] 0.3078365

6.2 Tunning do Modelo

Tunning ou Hypperparameter Optimization Este passo é importante para auxiliar a selecionar o melhor modelo. epsilon: 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 cost: 4, 8, 16, 32, 64, 128, 256, 512

ranges = list(epsilon = seq(0,1,.1), cost = 2^(2:9))
set.seed(123)
tmodel = tune(method = svm, WR_Grp ~ ., data = train2,
              ranges = ranges)
plot(tmodel)
summary(tmodel)

> summary(tmodel)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:
 epsilon cost
       0   16

- best performance: 0.2073372 

- Detailed performance results:
   epsilon cost     error dispersion
1      0.0    4 0.2159598 0.02054510
2      0.1    4 0.2159598 0.02054510
3      0.2    4 0.2159598 0.02054510
4      0.3    4 0.2159598 0.02054510
5      0.4    4 0.2159598 0.02054510
6      0.5    4 0.2159598 0.02054510
7      0.6    4 0.2159598 0.02054510
8      0.7    4 0.2159598 0.02054510
9      0.8    4 0.2159598 0.02054510
10     0.9    4 0.2159598 0.02054510
11     1.0    4 0.2159598 0.02054510
12     0.0    8 0.2114626 0.02161040
13     0.1    8 0.2114626 0.02161040
14     0.2    8 0.2114626 0.02161040
15     0.3    8 0.2114626 0.02161040
16     0.4    8 0.2114626 0.02161040
17     0.5    8 0.2114626 0.02161040
18     0.6    8 0.2114626 0.02161040
19     0.7    8 0.2114626 0.02161040
20     0.8    8 0.2114626 0.02161040
21     0.9    8 0.2114626 0.02161040
22     1.0    8 0.2114626 0.02161040
23     0.0   16 0.2073372 0.01888178
24     0.1   16 0.2073372 0.01888178
25     0.2   16 0.2073372 0.01888178
26     0.3   16 0.2073372 0.01888178
27     0.4   16 0.2073372 0.01888178
28     0.5   16 0.2073372 0.01888178
29     0.6   16 0.2073372 0.01888178
30     0.7   16 0.2073372 0.01888178
31     0.8   16 0.2073372 0.01888178
32     0.9   16 0.2073372 0.01888178
33     1.0   16 0.2073372 0.01888178
34     0.0   32 0.2095942 0.01475664
35     0.1   32 0.2095942 0.01475664
36     0.2   32 0.2095942 0.01475664
37     0.3   32 0.2095942 0.01475664
38     0.4   32 0.2095942 0.01475664
39     0.5   32 0.2095942 0.01475664
40     0.6   32 0.2095942 0.01475664
41     0.7   32 0.2095942 0.01475664
42     0.8   32 0.2095942 0.01475664
43     0.9   32 0.2095942 0.01475664
44     1.0   32 0.2095942 0.01475664
45     0.0   64 0.2125862 0.02271235
46     0.1   64 0.2125862 0.02271235
47     0.2   64 0.2125862 0.02271235
48     0.3   64 0.2125862 0.02271235
49     0.4   64 0.2125862 0.02271235
50     0.5   64 0.2125862 0.02271235
51     0.6   64 0.2125862 0.02271235
52     0.7   64 0.2125862 0.02271235
53     0.8   64 0.2125862 0.02271235
54     0.9   64 0.2125862 0.02271235
55     1.0   64 0.2125862 0.02271235
56     0.0  128 0.2148264 0.03043428
57     0.1  128 0.2148264 0.03043428
58     0.2  128 0.2148264 0.03043428
59     0.3  128 0.2148264 0.03043428
60     0.4  128 0.2148264 0.03043428
61     0.5  128 0.2148264 0.03043428
62     0.6  128 0.2148264 0.03043428
63     0.7  128 0.2148264 0.03043428
64     0.8  128 0.2148264 0.03043428
65     0.9  128 0.2148264 0.03043428
66     1.0  128 0.2148264 0.03043428
67     0.0  256 0.2253358 0.03007258
68     0.1  256 0.2253358 0.03007258
69     0.2  256 0.2253358 0.03007258
70     0.3  256 0.2253358 0.03007258
71     0.4  256 0.2253358 0.03007258
72     0.5  256 0.2253358 0.03007258
73     0.6  256 0.2253358 0.03007258
74     0.7  256 0.2253358 0.03007258
75     0.8  256 0.2253358 0.03007258
76     0.9  256 0.2253358 0.03007258
77     1.0  256 0.2253358 0.03007258
78     0.0  512 0.2298232 0.03273436
79     0.1  512 0.2298232 0.03273436
80     0.2  512 0.2298232 0.03273436
81     0.3  512 0.2298232 0.03273436
82     0.4  512 0.2298232 0.03273436
83     0.5  512 0.2298232 0.03273436
84     0.6  512 0.2298232 0.03273436
85     0.7  512 0.2298232 0.03273436
86     0.8  512 0.2298232 0.03273436
87     0.9  512 0.2298232 0.03273436
88     1.0  512 0.2298232 0.03273436

6.3 Best Model (Melhor Modelo)

Gráfico de performance do SVM. Onde a coloração em azul mais forte (escuro) aponta menor erro de classificação

mymodel = tmodel$best.model
summary(mymodel)

> summary(mymodel

Call:
best.tune(method = svm, train.x = WR_Grp ~ ., data = train2, ranges = ranges)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  16 
      gamma:  0.07692308 

Number of Support Vectors:  1308

 ( 646 662 )


Number of Classes:  2 

Levels: 
 [0, 6.5) 6.5, 10]

Matriz de Confusão e Erro de Classificação (Missclassification Error)

pred.bestsvm = predict(mymodel, train2)
(tab.bestsvm = table(Predito = pred.bestsvm, Real = train2$WR_Grp))
(erro.kernelRAD.BEST = 1 - sum(diag(tab.bestsvm))/sum(tab.bestsvm))

(tab.bestsvm = table(Predito = pred.bestsvm, Real = train2$WR_Grp))
          Real
Predito    [0, 6.5) 6.5, 10]
  [0, 6.5)     1172      236
  6.5, 10]      156     1103
> (erro.kernelRAD.BEST = 1 - sum(diag(tab.bestsvm))/sum(tab.bestsvm))
[1] 0.1469816

plot(mymodel, data = train2, num_voted_users~duration, 
     slice = list(title_year = 3, gross = 4))
plot(mymodel, data = train2, duration~num_voted_users, 
     slice = list(num_user_for_reviews = 3, gross = 4))

Gráf. de classificação do SVM

7 Conclusão e Considerações Finais

O nossos modelos, tanto RF quanto SVM Se ajustaram bem aos dados quanto à classificação filmes com alta e baixa aprovação.
O quadro a seguir apresenta o número de filmes em cada classificação segundo a base de teste, treino e inteira (movie). Aparentemente as categorias são bem balanceadas.

Base de dados	[0, 6.5)	[6.5, 10]
Movie (treino + teste)	1883	1883
Treino	1328	1339
Teste	555	544

Ainda conseguimos obter quais são os fatores/variáveis mais importantes para análise de classificação, dentre as variáveis existentes. Abaixo citamos as top 10 variáveis mais importantes.

hc6_1

Para o ajuste do RF

Utilizando a base de teste, o modelo retornou através da matriz de confusão uma sensibilidade de cerca de 85% e uma especificidade de 80%. Além disso, uma acurácia obtida de 83%.
É valido ressaltar, também, que o erro OOB estimado foi de 20%.
Ainda através da matriz de confusão, o nosso modelo erra em uma menor quantidade de vezes uma classificação de filmes acima de 6,5★ (baltaaprovação). A estimativa de erro de classificação desses filmes é de 16% versus 18% de erro de classificação em filmes com alta aprovação.
Abaixo podemos verificar a matriz de confusão do ajusto do RF (com base de teste), e afirmamos que o modelo erra menos na classificação dos filmes mais bem classificados:

PREDITO VS REAL	[0, 6.5)	[6.5, 10]	Erro de Classificação
[0, 6.5)	469	106	0.1843
[6.5, 10]	86	438	0.1641

Para o ajuste do SVM

O erro estimado foi de 14,698%, uma melhora de 5% do quando comparado com o ERRO OOB do RF.
Ainda através da matriz de confusão, o nosso modelo erra em uma menor quantidade de vezes uma classificação de filmes de acima 6,5★ (alta aprovação). A estimativa de erro de classificação desses filmes é de 12% versus 17% de erro de classificação em filmes com alta aprovação.
Abaixo podemos verificar a matriz de confusão do ajusto do RF (com base de teste), e afirmamos que o modelo erra menos na classificação dos filmes mais bem classificados:

PREDITO VS REAL	[0, 6.5)	[6.5, 10]	Erro de Classificação
[0, 6.5)	1172	236	0.1676
[6.5, 10]	156	1103	0.1239

Comparações

Erro de classficação	[0, 6.5)	[6.5, 10]
Random Forest (RF)	0.1843	0.1641
Support Vector Machine (SVM)	0.1676	0.1239

Em geral o SVM classificou melhor ambos os grupos de filmes (alto e baixo), com uma diferença de até 4% de acerto nas classificações em coparação com o RF na categoria de [6.5, 10].

Tempo de execução

proc.time() -ptm

##    user  system elapsed 
##   61.38    6.29   88.73

Modelo de classificação de notas de filmes

Uma aplicação com dados da plataforma IMDB - Internet Movie DataBases - via Random Forest (RF) e Support Vector Machine (SVM)

Ewerson C. Pimenta

Marcos Antônio E. de Oliveira

Rio de Janeiro, 29 de novembro de 2018